################################################################################
# Title:   instrument-preferences_data-analysis.R                                                    
# Script:  Data analysis for: "Understanding Policy Instrument Preferences under 
#          Conflicting Beliefs and Uncertainty"
# Date:    latest update on 25.04.2025                                                                             
# Authors: Milena Wiget (milena.wiget@eawag.ch), 
#.         Judit Lienert (judit.lienert@eawag.ch),
#          Karin Ingold (karin.ingold@unibe.ch)
# Project: TRAPEGO; WP 4.2 / 5.2
# Data:    (1) Survey "Akteursbefragung zu Entscheidungsprozessen zur Steuerung 
#              und Regulierung von Pflanzenschutzmitteln (2022)", conducted 
#              online in Limesurvey, May - Aug 2022
#          (2) Workshop on "Steuerung und Regulierung von Pflanzenschutzschutz-
#              mitteln in der Schweiz", focus group workshop at the University 
#              of Bern, 01.06.2023
# Input:    - instrument-preferences_model-data.csv
#           - instrument-preferences_weight-data.csv
# Output:   - model-results_BOLRMplotCust_C_wVar1.pdf --> Figures 2  
#           - model-results_BOLRMeffects_C_wVar1_R.pdf --> Figure 3 & 4
#           - data-results_cb-variance.pdf --> Figure A-1
#           - model-check_BOLRMfit_C_wVar1.pdf --> Figure A-2
#           - data-results_BOLRMest_C_wVar1_all.csv --> Table A-2
################################################################################

# R version 4.4.2 (2024-10-31) 
# Platform: aarch64-apple-darwin20 (64-bit)
# Running under: macOS Sequoia 15.4.1

# ==============================================================================
# Load libraries
# ==============================================================================

# list.functions.in.file("scripts/paper-preferences_data-analysis_B_toShare.R", alphabetic = TRUE)
lib_vect <- c("base","bayesplot","brms","corrplot","dplyr","forcats","ggeffects","ggplot2",
              "ggrepel","ggthemes","ggpubr","grDevices","grid","gridExtra","Hmisc",
              "marginaleffects","resample","rstan","rstatix","sjPlot","stats","tidybayes")
sapply(lib_vect,library,character.only=TRUE)

# ==============================================================================
# 1 Bayesian Ordinal Logistic Regression Model
# ==============================================================================
## 1.1 Data prep
## -----------------------------------------------------------------------------

# Load the data
BOLRMdata_basicVar <- as.data.frame(read.csv("instrument-preferences_model-data.csv", header=TRUE))
# Delete irrelevant columns
BOLRMdata_basicVarSel <- select(BOLRMdata_basicVar,-id)
# Set characters into factors
BOLRMdata_basicVarSel$policy <- factor(BOLRMdata_basicVarSel$policy,
                                    levels=c("Authorization","LimitValues",
                                    "UseRestriction","ZoneRestriction","Ban",
                                    "TaxRisk","TaxUpgrade","SubAppliTech",
                                    "SubFreeProduction","SubResearch",
                                    "InfoAccess","EarlyWarning","Training",
                                    "Agreements","Labels"))
BOLRMdata_basicVarSel$instrument <- factor(BOLRMdata_basicVarSel$instrument,
                                        levels=c("P","C","S","T","R"))
BOLRMdata_basicVarSel$type <- factor(BOLRMdata_basicVarSel$type,
                                  levels=c("AD","SC","IG"))


# Scale and round the continuous data
BOLRMdata_scaleVar <- BOLRMdata_basicVarSel
BOLRMdata_scaleVar[,c(4:12)] <- round(scale(BOLRMdata_basicVarSel[,c(4:12)]*100,scale=FALSE),1)
colnames(BOLRMdata_scaleVar) <- gsub("Var1","CB",as.character(colnames(BOLRMdata_scaleVar))) # CB for core belief

# Select data for the 5 policy instrument models 
BOLRMdata_Rvar <- dplyr::filter(BOLRMdata_scaleVar,instrument == "R")
BOLRMdata_Tvar <- dplyr::filter(BOLRMdata_scaleVar,instrument == "T")
BOLRMdata_Svar <- dplyr::filter(BOLRMdata_scaleVar,instrument == "S")
BOLRMdata_Cvar <- dplyr::filter(BOLRMdata_scaleVar,instrument == "C")
BOLRMdata_Pvar <- dplyr::filter(BOLRMdata_scaleVar,instrument == "P")

# Set dependent variable as factor variable for the ordinal linear regression model
BOLRMdata_Rvar$support <- factor(as.character(BOLRMdata_Rvar$support),
                              levels=c("1","2","3","4","5"),
                              labels=c("not at all","mostly not","neutral","mostly","fully"),
                              ordered = TRUE)
BOLRMdata_Rvar$policy <- factor(as.character(BOLRMdata_Rvar$policy),
                             levels=c("Authorization","LimitValues","UseRestriction",
                                      "ZoneRestriction","Ban"),
                             labels=c("A","LV","UR","ZR","B"))
BOLRMdata_Tvar$support <- factor(as.character(BOLRMdata_Tvar$support),
                              levels=c("1","2","3","4","5"),
                              labels=c("not at all","mostly not","neutral","mostly","fully"),
                              ordered = TRUE)
BOLRMdata_Tvar$policy <- factor(as.character(BOLRMdata_Tvar$policy),
                             levels=c("TaxRisk","TaxUpgrade"),
                             labels=c("TR","TU"))
BOLRMdata_Svar$support <- factor(as.character(BOLRMdata_Svar$support),
                              levels=c("2","3","4","5"),
                              labels=c("mostly not","neutral","mostly","fully"),
                              ordered = TRUE)
BOLRMdata_Svar$policy <- factor(as.character(BOLRMdata_Svar$policy),
                             levels=c("SubAppliTech","SubFreeProduction","SubResearch"),
                             labels=c("SAT","SFP","SR"))
BOLRMdata_Cvar$support <- factor(as.character(BOLRMdata_Cvar$support),
                              levels=c("2","3","4","5"),
                              labels=c("mostly not","neutral","mostly","fully"),
                              ordered = TRUE)
BOLRMdata_Cvar$policy <- factor(as.character(BOLRMdata_Cvar$policy),
                             levels=c("Agreements","Labels"),
                             labels=c("VA","L"))
BOLRMdata_Pvar$support <- factor(as.character(BOLRMdata_Pvar$support),
                              levels=c("3","4","5"),
                              labels=c("neutral","mostly","fully"),
                              ordered = TRUE)
BOLRMdata_Pvar$policy <- factor(as.character(BOLRMdata_Pvar$policy),
                             levels=c("InfoAccess","EarlyWarning","Training"),
                             labels=c("IA","EW","T"))

# List the data sets
BOLRMdataVar <- list(BOLRMdata_Rvar,BOLRMdata_Tvar,BOLRMdata_Svar,BOLRMdata_Cvar,BOLRMdata_Pvar)

## -----------------------------------------------------------------------------
## 1.2 Multilevel Ordinal Logistic Regression Model
## -----------------------------------------------------------------------------

# Estimate a Bayesian Multilevel Ordinal Logistic Regression model to specify 
# the relation between the dependent variable and explanatory variables)
# --> RA_Health is not included as it strongly correlates with RA_Environment
BOLRModel_C_wVar1 <- vector("list", length = 5)
names(BOLRModel_C_wVar1)[1:5] <- c("BOLRModel_C_wVar1_R","BOLRModel_C_wVar1_T",
                                   "BOLRModel_C_wVar1_S","BOLRModel_C_wVar1_C",
                                   "BOLRModel_C_wVar1_P")
for (i in 1:5){
  BOLRModel_C_wVar1[[i]] <- brm(formula = support ~ CB_Health + CB_Environment +
                                  CB_AgroEconomics + CB_SocioPolitics + 
                                  RA_Health + RA_Environment + 
                                  RA_AgroEconomics + RA_SocioPolitics + 
                                  CPSI + type + policy,
                              data = BOLRMdataVar[[i]], family= cumulative("logit"),
                              init = 0,iter=3000,seed = 123)
}

# Save results
BOLRMest_C_wVar1 <- vector("list", length = 5)
names(BOLRMest_C_wVar1)[1:5] <- c("BOLRMest_C_wVar1_R","BOLRMest_C_wVar1_T",
                                  "BOLRMest_C_wVar1_S","BOLRMest_C_wVar1_C",
                                  "BOLRMest_C_wVar1_P")
for (i in 1:5){
  BOLRMest_C_wVar1[[i]] <- summary(BOLRModel_C_wVar1[[i]])$fixed
  write.csv(BOLRMest_C_wVar1[[i]],paste0("data-results_",names(BOLRMest_C_wVar1)[i],".csv"),
            row.names = TRUE)
}


## -----------------------------------------------------------------------------
## 1.3 Model diagnostics
## -----------------------------------------------------------------------------

# Check the correlation among the variables
BOLRMcorm_C_wVar1 <- corrplot(cor(select(BOLRMdata_scaleVar,CB_Health,CB_Environment,
                                         CB_AgroEconomics,CB_SocioPolitics,
                                         RA_Health,RA_Environment, RA_AgroEconomics,
                                         RA_SocioPolitics)), 
                                  method = "ellipse", type = "lower", addCoef.col = 'black',
                                  tl.col = "black",tl.srt = 45)
pdf("model-check_BOLRMcorm_C_wVar1.pdf",18,18)
print(BOLRMcorm_C_wVar1)
dev.off()

# Check the chain plots
BOLRMchp_C_wVar1 <- vector("list", length = 5)
names(BOLRMchp_C_wVar1)[1:5] <- c("BOLRMchp_C_wVar1_R","BOLRMchp_C_wVar1_T",
                                  "BOLRMchp_C_wVar1_S","BOLRMchp_C_wVar1_C",
                                  "BOLRMchp_C_wVar1_P")
for (i in 1:5){
  BOLRMchp_C_wVar1[[i]] <- mcmc_trace(BOLRModel_C_wVar1[[i]])
  pdf(paste0("model-check_",names(BOLRMchp_C_wVar1)[i],".pdf"),18,18)
  print(mcmc_trace(BOLRModel_C_wVar1[[i]]))
  dev.off()
}

# Check the Rhat
BOLRMrhat_C_wVar1 <- vector("list", length = 5)
names(BOLRMrhat_C_wVar1)[1:5] <- c("BOLRMrhat_C_wVar1_R","BOLRMrhat_C_wVar1_T",
                                   "BOLRMrhat_C_wVar1_S","BOLRMrhat_C_wVar1_C",
                                   "BOLRMrhat_C_wVar1_P")
for (i in 1:5){
  BOLRMrhat_C_wVar1[[i]] <- mcmc_trace(BOLRModel_C_wVar1[[i]])
  pdf(paste0("model-check_",names(BOLRMrhat_C_wVar1)[i],".pdf"),24,12)
  print(mcmc_acf(BOLRModel_C_wVar1[[i]]))
  dev.off()
}

# Check the posterior predictive fit
BOLRMfit_C_wVar1 <- vector("list", length = 5)
names(BOLRMfit_C_wVar1)[1:5] <- c("BOLRMfit_C_wVar1_R","BOLRMfit_C_wVar1_T",
                                  "BOLRMfit_C_wVar1_S","BOLRMfit_C_wVar1_C",
                                  "BOLRMfit_C_wVar1_P")
for (i in 1:5){
  set.seed(123)
  y <- as.integer(BOLRMdataVar[[i]]$support)
  yrep <- posterior_predict(BOLRModel_C_wVar1[[i]])
  BOLRMfit_C_wVar1[[i]] <- ppc_bars(y,yrep)
}
# ... adjust graphical posterior predictive check (PPC) plots
BOLRMfit_C_wVar1[[1]] <- BOLRMfit_C_wVar1[[1]] + labs(x = "Support rating", fill="", color = "") + 
  scale_fill_manual(labels = "Observed data", values = "lightblue") +
  scale_colour_manual(labels = "Simulated data (posterior predictive distribution)", values = "darkslategrey") +
  theme(axis.title.x = element_text(face="bold",size=21, family = "sans",margin = margin(3,0,0,0,"mm")),
        axis.text.x = element_text(size= 21, family = "sans"),
        axis.title.y = element_text(face="bold",size=21, family = "sans",margin = margin(0,3,0,0,"mm")),
        axis.text.y = element_text(size= 21, family = "sans"),
        legend.position = c(0.35,0.9),
        legend.text = element_text(face="bold",size=21, family = "sans"))
BOLRMfit_C_wVar1[[2]] <- BOLRMfit_C_wVar1[[2]] + labs(x = "Support rating", fill="", color = "") + 
  scale_fill_manual(labels = "Observed data", values = "lightblue") +
  scale_colour_manual(labels = "Simulated data (posterior predictive distribution)", values = "darkslategrey") +
  theme(axis.title.x = element_text(face="bold",size=21, family = "sans",margin = margin(3,0,0,0,"mm")),
        axis.text.x = element_text(size= 21, family = "sans"),
        axis.title.y = element_text(face="bold",size=21, family = "sans",margin = margin(0,3,0,0,"mm")),
        axis.text.y = element_text(size= 21, family = "sans"),
        legend.position = c(0.35,0.9),
        legend.text = element_text(face="bold",size=21, family = "sans"))
BOLRMfit_C_wVar1[[3]] <- BOLRMfit_C_wVar1[[3]] + labs(x = "Support rating", fill="", color = "") + 
  expand_limits(x=0) +
  scale_x_continuous(breaks = c(0,1,2,3,4), labels = c("1","2","3","4","5")) +
  scale_fill_manual(labels = "Observed data", values = "lightblue") +
  scale_colour_manual(labels = "Simulated data (posterior predictive distribution)", values = "darkslategrey") +
  theme(axis.title.x = element_text(face="bold",size=21, family = "sans",margin = margin(3,0,0,0,"mm")),
        axis.text.x = element_text(size= 21, family = "sans"),
        axis.title.y = element_text(face="bold",size=21, family = "sans",margin = margin(0,3,0,0,"mm")),
        axis.text.y = element_text(size= 21, family = "sans"),
        legend.position = c(0.35,0.9),
        legend.text = element_text(face="bold",size=21, family = "sans"))
BOLRMfit_C_wVar1[[4]] <- BOLRMfit_C_wVar1[[4]] + labs(x = "Support rating", fill="", color = "") + 
  expand_limits(x=0) +
  scale_x_continuous(breaks = c(0,1,2,3,4), labels = c("1","2","3","4","5")) +
  scale_fill_manual(labels = "Observed data", values = "lightblue") +
  scale_colour_manual(labels = "Simulated data (posterior predictive distribution)", values = "darkslategrey") +
  theme(axis.title.x = element_text(face="bold",size=21, family = "sans",margin = margin(3,0,0,0,"mm")),
        axis.text.x = element_text(size= 21, family = "sans"),
        axis.title.y = element_text(face="bold",size=21, family = "sans",margin = margin(0,3,0,0,"mm")),
        axis.text.y = element_text(size= 21, family = "sans"),
        legend.position = c(0.35,0.9),
        legend.text = element_text(face="bold",size=21, family = "sans"))
BOLRMfit_C_wVar1[[5]] <- BOLRMfit_C_wVar1[[5]] + labs(x = "Support rating", fill="", color = "") + 
  expand_limits(x=c(-1,3)) +
  scale_x_continuous(breaks = c(-1,0,1,2,3), labels = c("1","2","3","4","5")) +
  scale_fill_manual(labels = "Observed data", values = "lightblue") +
  scale_colour_manual(labels = "Simulated data (posterior predictive distribution)", values = "darkslategrey") +
  theme(axis.title.x = element_text(face="bold",size=21, family = "sans",margin = margin(3,0,0,0,"mm")),
        axis.text.x = element_text(size= 21, family = "sans"),
        axis.title.y = element_text(face="bold",size=21, family = "sans",margin = margin(0,3,0,0,"mm")),
        axis.text.y = element_text(size= 21, family = "sans"),
        legend.position = c(0.35,0.9),
        legend.text = element_text(face="bold",size=21, family = "sans"))
# ... save PPC plots
for (i in 1:5){
  pdf(paste0("model-check_",names(BOLRMfit_C_wVar1)[i],".pdf"),12,6)
  print(BOLRMfit_C_wVar1[[i]])
  dev.off()
}

pdf("model-check_BOLRMfit_C_wVar1.pdf",24,18)
print(ggpubr::ggarrange(plotlist = BOLRMfit_C_wVar1,
                        labels = c("A","B","C","D","E"),
                        font.label = list(size = 21, face = "bold"),
                        ncol = 2, nrow = 3, align =c("hv"), heights = c(1,1,1,1,1)))
dev.off()

## -----------------------------------------------------------------------------
## 1.4 Result visualization
## -----------------------------------------------------------------------------
### 1.4.1 Model estimates
### ----------------------------------------------------------------------------

# Combine data frames with the results from the different models and save it
modelNames <- list("R","T","S","C","P")
categoryNames <- list("regulatory","market-based","market-based","cooperative","persuasive")
eVariables <- data.frame(var = c("CB_Health","CB_Environment","CB_AgroEconomics","CB_SocioPolitics",
                                   "RA_Health","RA_Environment","RA_AgroEconomics","RA_SocioPolitics"))
iVariables <- data.frame(var = c("Intercept[1]","Intercept[2]","Intercept[3]","Intercept[4]"))
c2Variables <- data.frame(var = c("policyLV","policyUR","policyZR","policyB",
                                   "policyTU","policySFP","policySR","policyL",
                                   "policyEW","policyT"))
c3Variables <- data.frame(var = c("typeSC","typeIG"))
DATAest_C_wVar1 <- BOLRMest_C_wVar1 
for (i in 1:5){
  DATAest_C_wVar1[[i]]$parameter <- rownames(BOLRMest_C_wVar1[[i]]) 
  DATAest_C_wVar1[[i]]$model <- as.vector(rep(modelNames[[i]],nrow(BOLRMest_C_wVar1[[i]])))
  DATAest_C_wVar1[[i]]$category <- as.vector(rep(categoryNames[[i]],nrow(BOLRMest_C_wVar1[[i]])))
  DATAest_C_wVar1[[i]]$parType <- ifelse(DATAest_C_wVar1[[i]]$parameter %in% eVariables$var, "explanatory",
                                         ifelse(DATAest_C_wVar1[[i]]$parameter %in% c2Variables$var, "controlP",
                                                ifelse(DATAest_C_wVar1[[i]]$parameter %in% iVariables$var, "intercept",
                                                       ifelse(DATAest_C_wVar1[[i]]$parameter %in% c3Variables$var, "controlA","controlC"))))
  rownames(DATAest_C_wVar1[[i]]) <- NULL
  DATAest_C_wVar1[[i]] <- select(DATAest_C_wVar1[[i]],category,model,parameter,parType,everything())
  DATAest_C_wVar1[[i]][,c(5:8)] <- round(DATAest_C_wVar1[[i]][,c(5:8)],2)
  DATAest_C_wVar1[[i]][,c(9:11)] <- round(DATAest_C_wVar1[[i]][,c(9:11)],0)
}
DATASETest_C_wVar1 <- bind_rows(DATAest_C_wVar1)
write.csv(DATASETest_C_wVar1,"data-results_BOLRMest_C_wVar1_all.csv", row.names = FALSE)

# Visualize the data sets from the different models 
# ... data preparation
DATASETvis_C_wVar1 <- DATASETest_C_wVar1 
DATASETvis_C_wVar1$parameter <- factor(DATASETvis_C_wVar1$parameter,
                                       levels=c("Intercept[1]","Intercept[2]","Intercept[3]","Intercept[4]",
                                               "CB_Health","CB_Environment","CB_AgroEconomics","CB_SocioPolitics",
                                               "RA_Health","RA_Environment","RA_AgroEconomics","RA_SocioPolitics",
                                               "CPSI","typeSC","typeIG","policyLV","policyUR","policyZR","policyB",
                                               "policyTU","policySFP","policySR","policyL","policyEW","policyT"),
                                       labels=c("Intercept 1","Intercept 2","Intercept 3","Intercept 4",
                                                "CB (Health)","CB (Env)","CB (Agro)","CB (Socio)",
                                                "RA (Health)","RA (Env)","RA (Agro)","RA (Socio)",
                                                "CPSI","Actor type (SC)","Actor type (IG)","Policy (LimitValues)",
                                                "Policy (useRestriction)","Policy (zoneRestriction)","Policy (Ban)",
                                                "Policy (TaxUpgrade)","Policy (DPfreeProduct)","Policy (SubResearch)",
                                                "Policy (Labels)","Policy (EarlyWarning)","Policy (Training)"))
DATASETvis_C_wVar1$parType <- factor(DATASETvis_C_wVar1$parType,
                                     levels=c("intercept","explanatory","controlC","controlA","controlP"),
                                     labels=c("Intercepts","Explanatory variables","Control I","Control II","Control III"))
DATASETvis_C_wVar1$category <- factor(DATASETvis_C_wVar1$category,
                                      levels=c("regulatory","market-based","cooperative","persuasive"))
DATASETvis_C_wVar1$model <- factor(DATASETvis_C_wVar1$model,
                                      levels=c("R","T","S","C","P"),
                                      labels=c("regulatory","market-based (tax)","market-based (subsidy)","cooperative","persuasive"))
colnames(DATASETvis_C_wVar1)[6:8] <- c("EstError","lowerCI95","upperCI95")
# ... plotting
pdf("model-results_BOLRMplotCust_C_wVar1.pdf",18,21)
plot <- ggplot(DATASETvis_C_wVar1, aes(x = Estimate, y = parameter)) + 
  geom_point(aes(color= model, shape = model), position=position_dodge(width=1), size=5) +
  geom_errorbar(aes(xmin=lowerCI95, xmax=upperCI95, color=model), 
                position=position_dodge(width=1), width=0, linewidth = 1.5) +
  geom_vline(xintercept=0,linetype='dashed', color='black', linewidth=0.5) + 
  theme_minimal() +
  labs(y="Model parameter",x="Model estimate (with 95% credibility interval)",color = "Instrument category:", shape = "Instrument category:") +
  scale_color_manual(name="Instrument category:", values = c("#500A60","#215274","#3E8B70","#D89C3C","#F5D35E")) +
  scale_shape_manual(name="Instrument category:", values = c(4,5,15,17,19)) +
  scale_y_discrete(limits = rev) +
  theme(axis.title.y = element_text(size = 21, face = "bold", margin=unit(c(0,0,0,0),"mm")),
        axis.text.y = element_text(size = 21),
        axis.title.x = element_text(size = 21, face ="bold",margin =margin(6,0,0,0,unit="mm")),
        axis.text.x.bottom = element_text(size = 21, margin =margin(2,0,0,0,unit="mm")),
        legend.title = element_text(size = 21, face = "bold",hjust=-1),
        legend.text = element_text(size = 21),
        legend.position = "top",
        legend.justification="right",
        legend.box.margin = margin(0,0,8,0,unit="mm")) +
  facet_wrap(parType~., ncol = 1, scales = "free",strip.position = "right") + 
  theme(strip.background = element_rect(color="white", fill="white", size=1.5, linetype="solid"),
        strip.text.y = element_text(size = 21, color = "black", face = "bold", margin =margin(6,0,0,3,unit="mm")),
        panel.spacing.y = unit(6, "mm"),
        strip.placement = "outside")
plotTable <- ggplotGrob(plot)
plotTable$heights[[9]] <- unit(4, "null")
plotTable$heights[[13]] <- unit(8, "null")
plotTable$heights[[17]] <- unit(1.6, "null")
plotTable$heights[[21]] <- unit(2, "null")
plotTable$heights[[25]] <- unit(3, "null")
grid.draw(plotTable)
dev.off()

### ----------------------------------------------------------------------------
### 1.4.2 Marginal Effects
### ----------------------------------------------------------------------------

# Create a list for the variable terms (continuous variables)
termsConV <- list("CB_Health","CB_Environment","CB_AgroEconomics","CB_SocioPolitics",
                  "RA_Health","RA_Environment","RA_AgroEconomics","RA_SocioPolitics",
                  "CPSI")

# Create a list for the x axis labels (continuous variables)
XaxisConV <- list("Relative weight of workload\n[deviation from the average]",
                  "Relative weight of NTO protection\n[deviation from the average]",
                  "Relative weight of economic viability\n[deviation from the average]",
                  "Relative weight of innovation potential\n[deviation from the average]",
                  "Risk attitude for health objectives\n[CE deviation from the average]",
                  "Risk attitude for environemental objectives\n[CE deviation from the average]",
                  "Risk attitude for agro-economic objectives\n[CE deviation from the average]",
                  "Risk attitude for socio-political objectives\n[CE deviation from the average]",
                  "CPSI [deviation from the average]")

# Create a level list for the different support levels
levelsAdj <- list(c("fully", "mostly", "neutral","mostly not","not at all"),
                  c("fully", "mostly", "neutral","mostly not","not at all"),
                  c("fully", "mostly", "neutral","mostly not"),
                  c("fully", "mostly", "neutral","mostly not"),
                  c("fully", "mostly", "neutral"))

# Create a title list for the models
titlesM <- list("Predictions of the marginal effect in the model for regulatory instruments",
                "Predictions of the marginal effect in the model for market-based tax incentives",
                "Predictions of the marginal effect in the model for market-based subsidies",
                "Predictions of the marginal effect in the model for cooperative instruments",
                "Predictions of the marginal effect in the model for persuasive instruments")

# Create a color list 
colorList <- list(c("#F5D35E","#D89C3C","#3E8B70","#215274","#500A60"),
                  c("#F5D35E","#D89C3C","#3E8B70","#215274","#500A60"),
                  c("#F5D35E","#D89C3C","#3E8B70","#215274"),
                  c("#F5D35E","#D89C3C","#3E8B70","#215274"),
                  c("#F5D35E","#D89C3C","#3E8B70"))

dashList <- list(c(1,5,2,4,3),c(1,5,2,4,3),c(1,5,2,4),c(1,5,2,4),c(1,5,2))
shapeList <- list(c(19,17,15,3,4),c(19,17,15,3,4),c(19,17,15,3),c(19,17,15,3),c(19,17,15))

# Create a list to save the marginal effect plots
BOLRMeffects_C_wVar1 <- vector("list", length = 5)
names(BOLRMeffects_C_wVar1)[1:5] <- c("BOLRMeffects_C_wVar1_R","BOLRMeffects_C_wVar1_T",
                                      "BOLRMeffects_C_wVar1_S","BOLRMeffects_C_wVar1_C",
                                      "BOLRMeffects_C_wVar1_P")

# Set the list names
for (i in 1:5) {
  BOLRMeffects_C_wVar1[[i]] <- vector("list", length = 10)
  names(BOLRMeffects_C_wVar1[[i]])[1:10] <- c(paste0(names(BOLRMeffects_C_wVar1)[i],"_CB_health"),
                                              paste0(names(BOLRMeffects_C_wVar1)[i],"_CB_env"),
                                              paste0(names(BOLRMeffects_C_wVar1)[i],"_CB_agro"),
                                              paste0(names(BOLRMeffects_C_wVar1)[i],"_CB_socio"),
                                              paste0(names(BOLRMeffects_C_wVar1)[i],"_RA_health"),
                                              paste0(names(BOLRMeffects_C_wVar1)[i],"_RA_env"),
                                              paste0(names(BOLRMeffects_C_wVar1)[i],"_RA_agro"),
                                              paste0(names(BOLRMeffects_C_wVar1)[i],"_RA_socio"),
                                              paste0(names(BOLRMeffects_C_wVar1)[i],"_CPSI"),
                                              paste0(names(BOLRMeffects_C_wVar1)[i],"_type"))
}

# Generate the marginal effect plots 
# ... for continuous variables
for (i in 1:5){
  for (j in 1:9){
    ggeff <- ggpredict(BOLRModel_C_wVar1[[i]], terms=termsConV[[j]]) 
    ggeff$response.level <- factor(ggeff$response.level,levels=levelsAdj[[i]])
    BOLRMeffects_C_wVar1[[i]][[j]] <- ggplot(ggeff, aes(x/100, predicted)) + 
      geom_line(aes(linetype=response.level, color=response.level), linewidth = 2) +
      geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=response.level), alpha=0.15) +
      theme_minimal() +
      labs(y="Predicted probability of support", x=XaxisConV[[j]],fill = "Support", 
           color = "Support", linetype = "Support") +
      scale_color_manual(values = colorList[[i]]) +
      scale_fill_manual(values = colorList[[i]]) +
      scale_linetype_manual(values = dashList[[i]]) +
      ylim(0.00,1.00) +
      theme(axis.title.x = element_text(face="bold",size=36, family = "sans",margin = margin(6,0,18,0,"mm")),
            axis.text.x = element_text(size= 36, family = "sans"),
            axis.title.y = element_text(face="bold",size=36, family = "sans",margin = margin(0,6,0,18,"mm")),
            axis.text.y = element_text(size= 36, family = "sans"),
            legend.title = element_text(size= 36, family = "sans", face = "bold"),
            legend.text = element_text(size= 36, family = "sans"),
            legend.position = c(0.2,0.8),
            legend.key.size = unit(1.5,"cm"),
            legend.background = element_rect(fill = "white",colour = 0))
  }
}
# ... for discrete variables
for (i in 1:5){
ggeff_t <- ggpredict(BOLRModel_C_wVar1[[i]], terms="type") 
ggeff_t$response.level <- factor(ggeff_t$response.level,levels=levelsAdj[[i]])
BOLRMeffects_C_wVar1[[i]][[10]] <- ggplot(ggeff_t, aes(x, predicted)) + 
  geom_point(aes(color=response.level, shape = response.level), position=position_dodge(width=0.4), size=7) +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high, color=response.level), 
                position=position_dodge(width=0.4), width=0.2, linewidth = 2) +
  theme_minimal() +
  labs(y="Predicted probability of support",x="Actor type",
       color = "Support", shape = "Support") +
  scale_color_manual(values = colorList[[i]]) +
  scale_shape_manual(values = shapeList[[i]]) +
  ylim(0.00,1.00) +
  theme(axis.title.x = element_text(face="bold",size=36, family = "sans",margin = margin(6,0,18,0,"mm")),
        axis.text.x = element_text(size= 36, family = "sans"),
        axis.title.y = element_text(face="bold",size=36, family = "sans",margin = margin(0,6,0,18,"mm")),
        axis.text.y = element_text(size= 36, family = "sans"),
        legend.title = element_text(size= 36, family = "sans", face = "bold"),
        legend.text = element_text(size= 36, family = "sans"),
        legend.position = c(0.95,0.8),
        legend.key.size = unit(1.5,"cm"),
        legend.background = element_rect(fill = "white",colour = 0))
}

# Save marginal effect plots 
for (i in 1:5){
  pdf(paste0("model-results_",names(BOLRMeffects_C_wVar1)[i],".pdf"),60,36)
  plot <- ggpubr::ggarrange(plotlist = BOLRMeffects_C_wVar1[[i]], 
                            labels = c("A","B","C","D","E","F","G","H","I","J"),
                            font.label = list(size = 36, face = "bold"),
                            ncol = 4, nrow = 3, align =c("hv"), heights = c(1,1,1))
  plotAn <- annotate_figure(plot, bottom = text_grob(titlesM[[i]], face = "bold", size = 36))
  print(plotAn)
  dev.off()
}

# Test correct visualization with plot_model function in R
# ... for the continuous variables
BOLRMeffects_test_conVar <- plot_model(BOLRModel_C_wVar1[[1]], type="pred",terms="CB_Environment") +
  labs(title = "", y="Predicted probability of support",x="Core Belief on the priority of environemental objectives (change in weight of the protection of NTO)") +
  theme(axis.title.x = element_text(face="bold",size=14, family = "sans",margin = margin(3,0,0,0,"mm")),
        axis.text.x = element_text(size= 12, family = "sans"),
        axis.title.y = element_text(face="bold",size=14, family = "sans",margin = margin(0,3,0,0,"mm")),
        axis.text.y = element_text(size= 12, family = "sans"))
BOLRMeffects_test_conVar$facet$params$nrow=1
BOLRMeffects_test_conVar
# ... for the discrete variables
BOLRMeffects_test_disVar <- plot_model(BOLRModel_C_wVar1[[1]], type="pred",terms="type") +
labs(title = "", y="Predicted probability of support",x="Core Belief on the priority of environemental objectives (change in weight of the protection of NTO)") +
  theme(axis.title.x = element_text(face="bold",size=14, family = "sans",margin = margin(3,0,0,0,"mm")),
        axis.text.x = element_text(size= 12, family = "sans"),
        axis.title.y = element_text(face="bold",size=14, family = "sans",margin = margin(0,3,0,0,"mm")),
        axis.text.y = element_text(size= 12, family = "sans"))
BOLRMeffects_test_disVar$facet$params$nrow=1
BOLRMeffects_test_disVar

# ==============================================================================
# 2 Variance
# ==============================================================================
## 2.1 Data preparation
## -----------------------------------------------------------------------------

# Load the data
matrix_wsGWsub <- as.matrix(select(as.data.frame(read.csv("instrument-preferences_weight-data.csv", header=TRUE)), -id))
colnames(matrix_wsGWsub) <- gsub("Q1_","",colnames(matrix_wsGWsub))
# Category files
health <- data.frame(obj = c("ApplicantsHealth","ConsumersHealth","Workload"))
environment <- data.frame(obj = c("NonTargetOrganisms","SoilProtection","ClimateProtection","WaterProtection"))
agroEconomics <- data.frame(obj = c("FarmersAutonomy","FoodSecurity","FoodCosts","EconomicViability"))
socioPolitics <- data.frame(obj = c("CostFairness","Coherence","InnovationPotential","JobSecurtiy","LandscapeQuality"))
# Calculate the variance of the relative importance of the policy core beliefs
variance_wsGWsub <- as.data.frame(colVars(matrix_wsGWsub,na.rm = TRUE))
colnames(variance_wsGWsub)[1] <- c("Var")
variance_wsGWsub$VarAd <- as.numeric(formatC(variance_wsGWsub$Var, format = "e", digits = 2))
# Add the sum of the relative weights of each policy core beliefs as indicator 
# of the overal importance
variance_wsGWsub$Sum <- colSums(matrix_wsGWsub,na.rm = TRUE)
# Add information on the domains of the policy core beliefs 
variance_wsGWsub$Domain <- ifelse(rownames(variance_wsGWsub) %in% health$obj, "Human health",
                                  ifelse(rownames(variance_wsGWsub) %in% environment$obj, "Environmental protection",
                                         ifelse(rownames(variance_wsGWsub) %in% agroEconomics$obj, "Agro-economic productivity",
                                                "Socio-political costs")))
# Add policy core belief labels
variance_wsGWsub$Labels <- rownames(variance_wsGWsub)


## -----------------------------------------------------------------------------
## 2.2 Visualization
## -----------------------------------------------------------------------------

# Visualize the weight variance and importance of the policy core beliefs across 
# all actors
pdf("data-results_cb-variance.pdf",18,9)
varScatterPlot <- ggplot(variance_wsGWsub, aes(x=VarAd, y=Sum)) + 
  geom_label_repel(aes(label=Labels, color=fct_inorder(Domain)), fontface = "bold", 
                   size =6, box.padding= 0.5, point.padding = 0.75,
                   segment.size= 1) +
  geom_point(aes(colour =fct_inorder(Domain)), size = 6) +
  theme_minimal()+
  labs(y="Overall importance (sum of the assigned relative weights)",x="Importance variance (variance of the assigned relative weights)",colour="Domain:") +
  scale_color_manual(values=c("#ADA7C9","#188C6F","#F5D35E","#005275")) +
  xlim(0.00e+00,1.5e-03) + ylim(0.0,2.5) +
  theme(axis.text.x=element_text(size=18,margin=margin(3,0,0,0,"mm")),
        axis.title.x=element_text(size=18,face="bold",margin=margin(6,0,0,0,"mm")),
        axis.text.y=element_text(size=18,margin=margin(0,3,0,0,"mm")),
        axis.title.y=element_text(size=18,face="bold",margin=margin(0,6,0,0,"mm")),
        legend.text=element_text(size=18),
        legend.title=element_text(size=18,face="bold"),
        legend.position = "top")
print(varScatterPlot)
dev.off()

## -----------------------------------------------------------------------------
# ==============================================================================


################################################################################ end of script